Hydrodynamic singularities and clustering in a freely cooling inelastic gas 
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We employ hydrodynamic equations to follow the clustering instability of a freely cooling dilute gas 
of inelastically colliding spheres into a well-developed nonlinear regime. We simplify the problem by 
dealing with a one-dimensional coarse-grained flow. We observe that at a late stage of the instability 
the shear stress becomes negligibly small, and the gas flows solely by inertia. As a result the flow 
formally develops a finite time singularity, as the velocity gradient and the gas density diverge 
at some location. We argue that flow by inertia represents a generic intermediate asymptotic of 
unstable free cooling of dilute inelastic gases. 

PACS numbers: 45.70.Qj, 47.70.Nd 
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A gas of inelastically colliding macroscopic particles is 
a simple paradigm of granular matter P, Q. IE 3 j and 
it appears in numerous applications, from astrophysics 
and geophysics to materials processing. One of the many 
fascinating phenomena in a freely cooling gas of inelas- 
tic particles is clustering instability |^, |a, S IM IM ^i 
im fl^ H^ . This instability attracted much attention 
from physicists to rapid granular flow 0,J3L, The clus- 
ters form an intricate cellular structure 0, 0, IE IE 0| ■ 
Though molecular dynamics (MD) simulations provide a 
valuable insight into the complicated dynamics of clus- 
tering, a better understanding requires a continuum the- 
ory. In this Letter we consider a freely cooling dilute gas 
of identical inelastic hard spheres. In this case a con- 
tinuum theory is derivable systematically from a kinetic 
equation, leading to hydrodynamic equations with a heat 
loss term caused by the inelasticity of particle collisions 
P, 0,0,0]. Hydrodynamics is expected to be an accu- 
rate leading order theory when the mean free path of the 
particles is much less than any length scale, and the mean 
time between two consecutive collisions is much less than 
any time scale, described hydrodynamically. Linearizing 
the hydrodynamic equations around a homogeneous cool- 
ing state (HCS), one finds two different linearly unstable 
modes: the shear mode and the thermal, or clustering, 
mode IE Si IE| ■ Growth of the shear mode corresponds to 
production of vorticity, while the clustering mode governs 
cluster formation. Nonlinear evolution of the clustering 
instability is a hard problem. Firstly, one has to deal 
here with nonlinear coupling of the shear and cluster- 
ing modes. Secondly, as the local density grows, the hy- 
drodynamic description becomes less accurate. It breaks 
down completely when the density approaches the point 
of the disorder-order transition in the gas of hard spheres. 

In this Letter we follow the clustering instability into 
a well-developed nonlinear stage by circumventing these 
two difhculties. Firstly, we put the particles into a long 
two-dimensional (2D) box, ^ Ly, so that the shear 
modes are strongly over-damped, and all coarse-grained 
quantities depend only on x. Secondly, we consider the 
limit of a very small area fraction of the particles. In this 
case, despite clustering, the gas density remains small. 



compared with the freezing density, for a very long time. 
Importantly, this limit does not preclude arbitrarily high 
density contrasts in the system. We solve the hydrody- 
namic equations numerically and observe that, at a late 
stage of the dynamics, the shear stress becomes negligi- 
bly small. As a result, the gas moves only by inertia, and 
the flow formally exhibits a finite time singularity. This 
singularity has a universal character if the initial mean 
velocity profile is smooth. We argue that flow by iner- 
tia is a generic intermediate asymptotic in more general 
multi-dimensional freely cooling granular flows, and that 
the finite time singularities form the skeleton of the later 
dynamics, when finite-density effects in the clusters come 
into play. 

Let each of N hard disks have a diameter a and mass 
TO = 1. Let the inelasticity of particle collisions be 
g = (1 — r)/2 > 0, where r is the (constant) coeffi- 
cient of normal restitution. Hydrodynamics deals with 
three coarse-grained fields: the number density n{x,t), 
the mean velocity v{x,t) and the granular tempera- 
ture T{x,t). We employ scaled variables n n/no, 



T -> T/To, 



1/2 



x/Lo- and t tT^^'^/L^, 



where no = N / {L^Ly) and Tq are the average number 
density and the initial temperature of the gas, respec- 
tively. In the dilute limit, granular hydrodynamic equa- 
tions 0,0,0113 read: 



dn dv ^ 
dt dx 



dv dP 

(a n— - — = , 
dt ox 



(b) (1) 



where d/dt ~ d/dt + vd/dx is the total time derivative, 
K — {2/y/n){aLxno)^^ is the Knudsen number which, 
up to a constant factor of order unity, is the ratio of the 
mean free path of the particles to L^, and P = —nT + 
{K/A) T^/'^{dv/dx) is the stress field. The validity of Eqs. 

and jSJ requires K <^ 1 (scale separation), ncr^ <C 1 
(dilute limit), and g <C 1 (nearly elastic collisions). 

The HCS is described by Half's cooling law T{x, t) — 
(1 + t/to)^'^, where to — K/Aq T|. A linear stability 
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analysis of the HCS, analogous to that of Refs. 
predicts clustering instability of the HCS if Kkx < 2g^/^, 
where is the (scaled) wave number of a small sinusoidal 
perturbation around the HCS. The wave number is quan- 
tized by the periodic boundary conditions: = 27rfc, 
where fc = 1, 2, . . . is the mode number. Therefore, the k- 
th mode is linearly unstable if irkK < q^l'^ The rest 
of the parameters fixed, the instability occurs when is 
sufficiently large. The number of the unstable modes in 
the system fcmax can serve as a measure of the instability 
magnitude. The growth/decay of small perturbations is 
algebraic. The density perturbations grow. The temper- 
ature and velocity perturbations decay, but the decays 
are slow compared to Haff 's cooling law. As a result, the 
flow tends to become supersonic 

We followed the clustering instability with fcmarr 3> 1 
into a strongly nonlinear regime by solving Eas. itlll and 
^ numerically. We used a Lagrangian scheme with 
periodic boundary conditions. The Lagrangian descrip- 
tion allowed us to resolve steep velocity gradients and 
high density peaks with good accuracy until close to sin- 
gularities, see below. The first series of hydrodynamic 
simulations dealt with generic initial conditions of the 
form n(x,t = 0) = H- 6n(x), T(x,t = 0) = 1 + bT{x) 
and v{x,t = 0) = Sv{x), where each of the small terms 
6n(x), 5T(x) and 5v{x) is a sum of a few hundred Fourier 
modes with random small amplitudes, of which a few 
dozen modes are linearly unstable. In all these simu- 
lations we observed strong clustering: development of 
multiple high and narrow density peaks, accompanied 
by steepening velocity gradients, as the gas temperature 
continues to decay. The gas density in the peaks grows 
without limit, until the time when our finite-difference 
scheme is unable to accurately follow the density growth 
in the highest density peak. The temporal growth of 
the density peaks, and of the velocity gradients, acceler- 
ates rapidly, implying a finite-time singularity. Figure 1 
shows a typical snapshot of the system close to the time 
of singularity. 

A convenient integral measure of the unstable cooling 
dynamics is the total energy of the system: 



E{t) 



1/2 / I 

[n T+—nv^ \ dx . 
-1/2 V 2 



(3) 



where nT is the thermal energy density, and nv'^/2 is 
the macroscopic kinetic energy density. A plot of E{t) is 
shown in Fig. 2. As expected, E{t) follows Haff's law at 
early times, but deviates from it at later times. Figure 2 
elucidates the role of each of the two terms in Eq. 
Both the thermal energy, and the macroscopic kinetic en- 
ergy initially decay with time; the thermal energy decays 
faster. At later times the kinetic energy approaches a 
constant. As a result, E{t) is dominated by the ther- 
mal energy at early times and by the kinetic energy at 
later times. Remarkably, the thermal energy continues 
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FIG. 1: The density and velocity profiles at scaled time 
t = 7.043, shortly before the major density peak develops 
singularity. The parameters K = A - 10"'' and q = 10~^ corre- 
spond to 79 linearly unstable Fourier modes. 2000 Lagrangian 
mesh points are used, so the major density peak includes more 
than 50 mesh points above the density value of n = 10^ . The 
inset shows an earlier density history (at t = 2, 3 and 4) of a 
region around the major density peak. 



to follow Haff's law until the time of singularity. 




FIG. 2: The total energy of the system E (thick solid line), 
the thermal energy (circles), the macroscopic kinetic energy 
(squares) and Haff's cooling law (thin solid line) versus time 
for the simulation shown in Fig. 1. 

The finite-time singularities of the velocity gradient 
and the density strongly indicate that, at late times, the 
stress tensor P becomes negligibly small, and the gas 
flows by inertia only. An additional evidence for the flow 
by inertia is the constancy of the macroscopic kinetic 
energy at late times. The flow by inertia is described by 
the equation 

dv dv ^ 
dt'-''d-x=' 
and Eq. JQla. This problem is soluble analytically [20| : 



v{x,t) ^ Vq{^) , (a) n{x,t) 



(b) (5) 



where Vq{£_) = dvniO/dS., while wo(^) and no{^) are the 
velocity and density of the gas, respectively, at some "ini- 
tial" moment of time (which should be late enough so 



that the flow by inertia has already set in). The relation 
between Eulerian coordinate x and Lagrangian coordi- 
nate ^ is the following: x — + vq{^) t. The finite-time 
singularities of both the velocity gradient 



dv{x, t) 
dx 



(6) 



and the density, Eq. ©b, occur when the denominator 
in Eq. © becomes zero for the first time. We com- 
pared these predictions with a numerical solution of the 
full hydrodynamic equations (Q) and (0), for the same 
parameters K = 4- 10""* and q ~ 10^^, but with simpler 
initial conditions: n{x,t = 0) = T{x,t = 0) = 1, and a 
single Fourier mode for the velocity: 



v{x, t = 0) ^ a sin(27rx) 



-0.05. 



(7) 



In this case only one singularity develops (at x — 0). 
Figure 3 shows the gas velocity v versus ^ = x — tv{x,t) 
at different times. The different curves collapse into a 
single curve with an accuracy better than 1.5%. Ad- 
ditional tests deal with the behavior of the flow in the 
close vicinity of x = 0, as the singularity time is ap- 
proached. For this smooth symmetric flow we can write 
voiO = -^/t + + where t = r is the time 

of singularity in the flow-by-inertia model, and C > 
is a constant. In the Eulerian coordinates this yields a 
solution in an implicit form. In the leading order 



-t'v{x,t')-CT^v^{x,t') 



(8) 



where t' = r — i is the time to the singularity. Not 
too close to the singularity point a; = one obtains 

V ^ {—xY^^. As the velocity profile lO is self-similar: 
v{x,t') it')'^/'^V [x/{t')^/'^], the velocity gradient is 
dv/dx ~ {t')^-'-dV/dw, where w — x/{t')^^'^. The shape 
function is determined by the equation Ct'^V^ + 

V + w = 0. What is the density behavior close to the 
singularity? Very close to a; = the density grows in- 
definitely: no(0)(l — t/r)^^; outside of that region (but 
still close enough to x — Q) n{x, t) approaches a univer- 



sal profile n 
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,21] . We verified these properties 



numerically, see examples in Fig. 4. Importantly, for a 
strong instability, kmax ^ 1, the system "freezes up", 
and the motion by inertia sets in very rapidly. Indeed, 
the scaled velocity profile in Fig. 3 is very close to the 
initial profile Q . Building on this simplification, we can 
expand Eq. Q in the vicinity of a: = and predict the 
time of singularity: r = (27r|a|)^^ ~ 3.18 which agrees 
within 2% with the simulation result, see Fig. 4a. In 
addition, the linear time dependence of the quantities, 
shown in the inset of Fig. 4a, sets in already at early 
times. 

Therefore, a strongly nonlinear regime of the quasi- 
one-dimensional clustering instability in a dilute granular 
flow is describable by a simple flow by inertia, until the 
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FIG. 3: The numerically computed velocity is shown versus 
X (a) and versus ^ = x — vt (b) at times 1, 2 and 3.225 
(the profiles in figure a steepen as the time progresses). Also 
shown in figure b is the initial profile Q. All the curves in 
figure b coincide within 1.5%. The simulation parameters are 
K = 4-10-'^ andg = 10-^ 



200t 
£150 

X 

5 100 
50 

c 

0* 



_ 4 
« 3 
I 2 




(a) 



12 3 4 



10^ 
10^ 
10^ 
10° 



(b) 



io''io''ioMo' 



FIG. 4: Numerically computed values of \dv/dx\ (filled 
squares) and n (empty circles) at x=0 versus time (a). The 
inset shows the respective inverse values. Figure b depicts the 
spatial profiles of n at time moments 3.209 and 3.225. The 
straight line is a x~'^^^ dependence; it is given for reference. 
The parameters are the same as in Fig. 3. 



moment of singularity |22|. In a related work Ben-Naim 
et al. jl^ investigated the dynamics of point-like par- 
ticles, inelastically colliding on a line. The strictly one- 
dimensional setting of Ref. |23| makes a hydrodynamic 
description problematic. Still, Ben-Naim et al. observed 
that the Burgers equation with vanishing viscosity is a 
proper continuum model for their system. It remains to 
be seen whether the Burgers equation or some other sat- 
uration mechanism applies to our gitosi-one-dimensional 
model at a later stage of the dynamics, when finite- 
density effects come into play. We stress that the (hydro- 
dynamic) density singularities are entirely different from 
inelastic collapse l24|| (divergence of the particle collision 
rate at some locations) which is a discrete-particle effect. 
We also note in passing that statistical properties of the 
flow by inertia (for example, the dynamics of the struc- 
ture function) are well understood [25|. 

What can be said about a fully multi-dimensional 
strongly unstable cooling flow, when unstable shear and 
clustering modes are coupled? A natural assumption, 
motivated already by the linear theory of the cluster- 
ing/shearing instabilities of the HCS 0, IIIj is that 
the stress tensor "freezes up" , and flow by inertia sets in 
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here as well. A possible counter- argument involves vis- 
cous heating of the system by the unstable shear modes. 
The heating effect is absent in the linear regime of the 
instability (as the viscous heating is of the second order 
with respect to the perturbation amplitude), but it comes 
into play in the nonlinear regime. The present state of 
theory makes it difficult to prove that the viscous heating 
cannot arrest, in some locations, the freezing of the stress 
tensor. However, MD simulations in 2D strongly indicate 
that the freezing continues unarrested. For example, Nie 
et al. 01 observe that, at late times, "the thermal energy 
becomes much smaller than the (macroscopic) kinetic en- 
ergy" . Based on this evidence we argue that, in the dilute 
regime, this strongly supersonic flow should be describ- 
able by muUi- dimensional flow-by-inertia equations 

av/9t + (v • V)v = , 9n/at + V • (nv) = 0. (9) 

This flow also exhibits finite-time singularities |2ll Efij i , 
and the singularities form cellular structures, most of 
the material being concentrated along the cell bound- 
aries This picture resembles the density distribu- 
tion of granular clusters observed in MD simulations 
of freely cooling gases of inelastic hard spheres in 2D 
Interestingly, the multi-dimensional singu- 
larities of Eqs. were studied previously in an entirely 
different context: in the so called Zeldovich approxima- 
tion of theory of formation of structure in an expanding 
universe [2H . 

We stress that there are important differences between 
the multi-dimensional clustering instability and the Zel- 
dovich model. In the process of clustering instability of 
inelastic gases a considerable vorticity is generated, while 
in Zeldovich model the flow is assumed to be potential 
|2ll |. Still, it was found, in a rare treatment of a more 
general (non-potential) velocity field, that "high-density 
regions should be high-vorticity regions' ' This find- 
ing appears to agree with MD simulations of freely cool- 
ing granular gases in 2D 

In summary, by following the unstable cooling dynam- 
ics of a dilute inelastic gas we identified an important new 
intermediate asymptotic regime: a nonlinear flow by iner- 
tia. We argue that high-density regions in the gas, which 
are precursors of densely packed granular clusters, are 
caused by the flow by inertia, rather than directly by the 
pressure gradient. Our results indicate that the role of 
the clustering and shearing instabilities of the free cooling 
is "merely" to produce a long-lived spatially non-uniform 
supersonic velocity field needed for the development of 
the high-density regions by the flow by inertia. There- 
fore, a due account of the flow-by-inertia regime will be 
important in the future theory of "life after singularity" , 
where the singularities are smoothed by finite-density ef- 
fects in the clusters, and a coarsening process develops 
. No first-principles coarse-grained description of that 
final stage is yet available. 
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